if(t > 1) y[t] <- alpha*y[t-1] + beta*x[t] + u[t]
}
###Prepare time series and lagged variable data
ynl <- y[-(1:2)]
ylag1 <- y[-c(1, length(y))]
ylag2 <- y[-(length(y):(length(y)-1)) ]
xn <- x[-1]
xlag <- x[-length(x)]
xn <- xn[-1]
xlag <- xlag[-1]
eq4 <- lm(ynl ~ 0 + ylag1 + ylag2 + xn + xlag)
lgdv <- lm(ynl ~ 0 + ylag1 + xn)
lgdv2 <- lm(ynl ~ 0 + ylag1 + ylag2 + xn)
reg <- lm(ynl ~ 0 + xn)
# EQ4
mat_coef[[j]][i, 1] <- summary(eq4)$coefficients[3, 1]
mat_se[[j]][i, 1] <- summary(eq4)$coefficients[3, 2]
mat_pvals[[j]][i, 1] <- summary(eq4)$coefficients[3, 4]
mat_lb[[j]][i, 1] <- confint(eq4)[3, 1]
mat_ub[[j]][i, 1] <- confint(eq4)[3, 2]
# LGDV
mat_coef[[j]][i, 2] <- summary(lgdv)$coefficients[2, 1]
mat_se[[j]][i, 2] <- summary(lgdv)$coefficients[2, 2]
mat_pvals[[j]][i, 2] <- summary(lgdv)$coefficients[2, 4]
mat_lb[[j]][i, 2] <- confint(lgdv)[2, 1]
mat_ub[[j]][i, 2] <- confint(lgdv)[2, 2]
# LGDV2
mat_coef[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 1]
mat_se[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 2]
mat_pvals[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 4]
mat_lb[[j]][i, 3] <- confint(lgdv2)[3, 1]
mat_ub[[j]][i, 3] <- confint(lgdv2)[3, 2]
# REG
mat_coef[[j]][i, 4] <- summary(reg)$coefficients[1, 1]
mat_se[[j]][i, 4] <- summary(reg)$coefficients[1, 2]
mat_pvals[[j]][i, 4] <- summary(reg)$coefficients[1, 4]
mat_lb[[j]][i, 4] <- confint(reg)[1, 1]
mat_ub[[j]][i, 4] <- confint(reg)[1, 2]
}
}
View(mat_coef)
View(eq4)
View(mat_coef)
mat_coef[[1]][1]
mat_coef[[1]][, 1]
## phi2 = 0
### EQ4
plot(density(mat_coef[[1]][, 1]))
shapiro.test(mat_coef[[1]][, 1])
### LGDV
plot(density(mat_coef[[1]][, 2]))
shapiro.test(mat_coef[[1]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[1]][, 3]))
shapiro.test(mat_coef[[1]][, 3]) # not normal
### REG
plot(density(mat_coef[[1]][, 4]))
shapiro.test(mat_coef[[1]][, 4]) # not normal
## EQ4, phi2 = 0.1
### EQ4
plot(density(mat_coef[[2]][, 1]))
shapiro.test(mat_coef[[2]][, 1]) # not normal
### LGDV
plot(density(mat_coef[[2]][, 2]))
shapiro.test(mat_coef[[2]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[2]][, 3]))
shapiro.test(mat_coef[[2]][, 3]) # not normal
### REG
plot(density(mat_coef[[2]][, 4]))
shapiro.test(mat_coef[[2]][, 4]) # not normal
## EQ4, phi2 = 0.2
### EQ4
plot(density(mat_coef[[3]][, 1]))
shapiro.test(mat_coef[[3]][, 1]) # not normal at 90%
### LGDV
plot(density(mat_coef[[3]][, 2]))
shapiro.test(mat_coef[[3]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[3]][, 3]))
shapiro.test(mat_coef[[3]][, 3]) # not normal
### REG
plot(density(mat_coef[[3]][, 4]))
shapiro.test(mat_coef[[3]][, 4]) # not normal
## EQ4, phi2 = 0.3
### EQ4
plot(density(mat_coef[[4]][, 1]))
shapiro.test(mat_coef[[3]][, 1]) # normal
shapiro.test(mat_coef[[4]][, 1]) # normal
### LGDV
plot(density(mat_coef[[4]][, 2]))
shapiro.test(mat_coef[[4]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[4]][, 3]))
shapiro.test(mat_coef[[4]][, 3]) # not normal
### REG
plot(density(mat_coef[[4]][, 4]))
shapiro.test(mat_coef[[4]][, 4]) # not normal
## EQ4, phi2 = 0.4
### EQ4
plot(density(mat_coef[[5]][, 1]))
shapiro.test(mat_coef[[5]][, 1]) # normal
### LGDV
plot(density(mat_coef[[5]][, 2]))
shapiro.test(mat_coef[[5]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[5]][, 3]))
shapiro.test(mat_coef[[5]][, 3]) # not normal
### REG
plot(density(mat_coef[[5]][, 4]))
shapiro.test(mat_coef[[5]][, 4]) # not normal
# Perecent Bias
bias <- matrix(NA, 6, 4)
## EQ4, phi2 = 0
### EQ4
plot(density(mat_coef[[6]][, 1]))
shapiro.test(mat_coef[[6]][, 1]) # normal
### LGDV
plot(density(mat_coef[[6]][, 2]))
shapiro.test(mat_coef[[6]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[6]][, 3]))
shapiro.test(mat_coef[[6]][, 3]) # not normal
### REG
plot(density(mat_coef[[6]][, 4]))
shapiro.test(mat_coef[[6]][, 4]) # not normal
load("/Users/alikagalwala/Dropbox/CSTS-Panel Papers/JOP RRR/Replication Files/Clark and Linzer 2015/Sluggish/2000sims_sluggish.RData")
normal.pool <- matrix(data = NA, nrow = 99, ncol = 6)
colnames(normal.pool) <- c("units", "within_units", "corr", "p_value", "normal_90", "normal_95")
normal.fe <- matrix(data = NA, nrow = 99, ncol = 6)
colnames(normal.fe) <- c("units", "within_units", "corr", "p_value", "normal_90", "normal_95")
normal.re <- matrix(data = NA, nrow = 99, ncol = 6)
colnames(normal.re) <- c("units", "within_units", "corr", "p_value", "normal_90", "normal_95")
index <- 1
j <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
for (i in 1:length(truecor)) {
# pooled
normal.pool[index, 1] <- unitloop
normal.pool[index, 2] <- Nperloop
normal.pool[index, 3] <- truecor[i]
normal.pool[index, 4] <- shapiro.test(coeflist.pool[[j]][ , i])$p.value
normal.pool[index, 5] <- ifelse(shapiro.test(coeflist.pool[[j]][ , i])$p.value > 0.1, "normal", "not normal")
normal.pool[index, 6] <- ifelse(shapiro.test(coeflist.pool[[j]][ , i])$p.value > 0.05, "normal", "not normal")
# fe
normal.fe[index, 1] <- unitloop
normal.fe[index, 2] <- Nperloop
normal.fe[index, 3] <- truecor[i]
normal.fe[index, 4] <- shapiro.test(coeflist.fe[[j]][ , i])$p.value
normal.fe[index, 5] <- ifelse(shapiro.test(coeflist.fe[[j]][ , i])$p.value > 0.1, "normal", "not normal")
normal.fe[index, 6] <- ifelse(shapiro.test(coeflist.fe[[j]][ , i])$p.value > 0.05, "normal", "not normal")
# re
normal.re[index, 1] <- unitloop
normal.re[index, 2] <- Nperloop
normal.re[index, 3] <- truecor[i]
normal.re[index, 4] <- shapiro.test(coeflist.re[[j]][ , i])$p.value
normal.re[index, 5] <- ifelse(shapiro.test(coeflist.re[[j]][ , i])$p.value > 0.1, "normal", "not normal")
normal.re[index, 6] <- ifelse(shapiro.test(coeflist.re[[j]][ , i])$p.value > 0.05, "normal", "not normal")
index <- index + 1
print(j)
}
j <- j + 1
}
}
View(normal.fe)
View(normal.pool)
View(normal.re)
##### ---- Figure 1c and 1d Replication ---- #####
rm(list = ls())
# Load required libraries
library(tidyverse)
library(cowplot)
library(robustbase)
library(matrixStats)
# Set working directory
#setwd("/Users/alikagalwala/Dropbox/CSTS-Panel Papers/JOP RR/Replication Files/Wilkins 2018")
setwd("C:/Users/alikagalwala/Dropbox/CSTS-Panel Papers/JOP RRR/Replication Files/Wilkins 2018")
# Set seed
set.seed(1381) # Same as Wilkins
# Parameters for the DGP
phi2 <- seq(0, 0.5, 0.1)
beta <- 0.5
alpha <- 0.75
rho <- 0.95
n <- 100
sims <- 1000 # No. of MC sims
# List to store results
mat_coef <- list()
mat_pvals <- list()
mat_lb <- list()
mat_ub <- list()
mat_se <- list()
# Monte carlo analysis, figure 1c and 1d
for(j in 1:length(phi2)){
phi <- phi2[j]
# Matrices to store results
mat_coef[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
colnames(mat_coef[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
mat_se[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
colnames(mat_se[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
mat_pvals[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
colnames(mat_pvals[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
mat_lb[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
colnames(mat_lb[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
mat_ub[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
colnames(mat_ub[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
for(i in 1:sims) {
# Generate autoregressive x
x <- numeric(n+2)
x[1] <- rnorm(1, 0, 1)
for(t in 2:length(x)){
x[t] <- rho*x[t-1] + rnorm(1, 0, 1)
}
# Generate u (serially correlated errors)
u <- numeric(n+2)
u[1] <- rnorm(1, 0, 1)
for(t in 2:length(u)){
u[t] <- phi*u[t-1] + rnorm(1, 0, 1)
}
# Generate y
y <- numeric(n+2)
yinit <- rnorm(1, 0, 1)
for(t in 1:length(x)){
if(t==1) y[t] <- alpha*yinit + beta*x[t] + u[t]
if(t > 1) y[t] <- alpha*y[t-1] + beta*x[t] + u[t]
}
###Prepare time series and lagged variable data
ynl <- y[-(1:2)]
ylag1 <- y[-c(1, length(y))]
ylag2 <- y[-(length(y):(length(y)-1)) ]
xn <- x[-1]
xlag <- x[-length(x)]
xn <- xn[-1]
xlag <- xlag[-1]
eq4 <- lm(ynl ~ 0 + ylag1 + ylag2 + xn + xlag)
lgdv <- lm(ynl ~ 0 + ylag1 + xn)
lgdv2 <- lm(ynl ~ 0 + ylag1 + ylag2 + xn)
reg <- lm(ynl ~ 0 + xn)
# EQ4
mat_coef[[j]][i, 1] <- summary(eq4)$coefficients[3, 1]
mat_se[[j]][i, 1] <- summary(eq4)$coefficients[3, 2]
mat_pvals[[j]][i, 1] <- summary(eq4)$coefficients[3, 4]
mat_lb[[j]][i, 1] <- confint(eq4)[3, 1]
mat_ub[[j]][i, 1] <- confint(eq4)[3, 2]
# LGDV
mat_coef[[j]][i, 2] <- summary(lgdv)$coefficients[2, 1]
mat_se[[j]][i, 2] <- summary(lgdv)$coefficients[2, 2]
mat_pvals[[j]][i, 2] <- summary(lgdv)$coefficients[2, 4]
mat_lb[[j]][i, 2] <- confint(lgdv)[2, 1]
mat_ub[[j]][i, 2] <- confint(lgdv)[2, 2]
# LGDV2
mat_coef[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 1]
mat_se[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 2]
mat_pvals[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 4]
mat_lb[[j]][i, 3] <- confint(lgdv2)[3, 1]
mat_ub[[j]][i, 3] <- confint(lgdv2)[3, 2]
# REG
mat_coef[[j]][i, 4] <- summary(reg)$coefficients[1, 1]
mat_se[[j]][i, 4] <- summary(reg)$coefficients[1, 2]
mat_pvals[[j]][i, 4] <- summary(reg)$coefficients[1, 4]
mat_lb[[j]][i, 4] <- confint(reg)[1, 1]
mat_ub[[j]][i, 4] <- confint(reg)[1, 2]
}
}
# Normality tests
## phi2 = 0
### EQ4
plot(density(mat_coef[[1]][, 1]))
shapiro.test(mat_coef[[1]][, 1]) # not normal
### LGDV
plot(density(mat_coef[[1]][, 2]))
shapiro.test(mat_coef[[1]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[1]][, 3]))
shapiro.test(mat_coef[[1]][, 3]) # not normal
### REG
plot(density(mat_coef[[1]][, 4]))
shapiro.test(mat_coef[[1]][, 4]) # not normal
## EQ4, phi2 = 0.1
### EQ4
plot(density(mat_coef[[2]][, 1]))
shapiro.test(mat_coef[[2]][, 1]) # not normal at 90%
### LGDV
plot(density(mat_coef[[2]][, 2]))
shapiro.test(mat_coef[[2]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[2]][, 3]))
shapiro.test(mat_coef[[2]][, 3]) # not normal
### REG
plot(density(mat_coef[[2]][, 4]))
shapiro.test(mat_coef[[2]][, 4]) # not normal
## EQ4, phi2 = 0.2
### EQ4
plot(density(mat_coef[[3]][, 1]))
shapiro.test(mat_coef[[3]][, 1]) # normal
### LGDV
plot(density(mat_coef[[3]][, 2]))
shapiro.test(mat_coef[[3]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[3]][, 3]))
shapiro.test(mat_coef[[3]][, 3]) # not normal
### REG
plot(density(mat_coef[[3]][, 4]))
shapiro.test(mat_coef[[3]][, 4]) # not normal
## EQ4, phi2 = 0.3
### EQ4
plot(density(mat_coef[[4]][, 1]))
shapiro.test(mat_coef[[4]][, 1]) # normal
### LGDV
plot(density(mat_coef[[4]][, 2]))
shapiro.test(mat_coef[[4]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[4]][, 3]))
shapiro.test(mat_coef[[4]][, 3]) # not normal
### REG
plot(density(mat_coef[[4]][, 4]))
shapiro.test(mat_coef[[4]][, 4]) # not normal
## EQ4, phi2 = 0.4
### EQ4
plot(density(mat_coef[[5]][, 1]))
shapiro.test(mat_coef[[5]][, 1]) # normal
### LGDV
plot(density(mat_coef[[5]][, 2]))
shapiro.test(mat_coef[[5]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[5]][, 3]))
shapiro.test(mat_coef[[5]][, 3]) # not normal
### REG
plot(density(mat_coef[[5]][, 4]))
shapiro.test(mat_coef[[5]][, 4]) # not normal
## EQ4, phi2 = 0.5
### EQ4
plot(density(mat_coef[[6]][, 1]))
shapiro.test(mat_coef[[6]][, 1]) # normal
### LGDV
plot(density(mat_coef[[6]][, 2]))
shapiro.test(mat_coef[[6]][, 2]) # not normal
### LGDV2
plot(density(mat_coef[[6]][, 3]))
shapiro.test(mat_coef[[6]][, 3]) # not normal
### REG
plot(density(mat_coef[[6]][, 4]))
shapiro.test(mat_coef[[6]][, 4]) # not normal
# Perecent Bias
bias <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
bias[i, ] <- (colMeans(mat_coef[[i]] - beta))/beta * 100
}
bias <- cbind(phi2, bias)
colnames(bias)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
bias <- as.tibble(bias)
bias <- bias %>%
select(-REG) %>%
gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "bias")
percentbias <- ggplot(bias, aes(x = phi, y = bias, color = Model)) +
geom_line(aes(linetype = Model), size = 1) +
ylim(-50, 50) +
ylab("") +
xlab(expression(phi)) +
ggtitle("Percent Bias") +
theme_minimal() +
theme(legend.position = "none") +
scale_linetype_manual(values=c("solid", "longdash", "dotdash"))
# RMSE
rmse <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
rmse[i, ] <- sqrt(colMeans((mat_coef[[i]] - beta)^2))
}
rmse <- cbind(phi2, rmse)
colnames(rmse)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
rmse <- as.tibble(rmse)
rmse <- rmse %>%
select(-REG) %>%
gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "rmse")
rootmean <- ggplot(rmse, aes(x = phi, y = rmse, color = Model)) +
geom_line(aes(linetype = Model), size = 1) +
ylab("") +
ylim(0, 0.25) +
xlab(expression(phi)) +
ggtitle("RMSE") +
theme_minimal() +
theme(legend.position = "none") +
scale_linetype_manual(values=c("solid", "longdash", "dotdash"))
# # Mean Absolute Bias
# mab <- matrix(NA, 6, 4)
# for(i in 1:length(phi2)){
#   mab[i, ] <- colMeans(abs(mat_coef[[i]] - beta))
# }
# mab <- cbind(phi2, mab)
# colnames(mab)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
# mab <- as.tibble(mab)
# mab <- mab %>%
#   select(-REG) %>%
#   gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "mab")
#
# meanabs <- ggplot(mab, aes(x = phi, y = mab, color = Model)) +
#   geom_line(aes(linetype = Model)) +
#   ylim(0, 0.2) +
#   ylab("") +
#   xlab(expression(phi)) +
#   ggtitle("Mean Absolute Bias") +
#   theme(legend.position = "none")
# Power
power <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
power[i, ] <- colMeans(ifelse(mat_pvals[[i]] <= 0.05, 1, 0))
}
power <- cbind(phi2, power)
colnames(power)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
power <- as.tibble(power)
power <- power %>%
select(-REG) %>%
gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "power")
pow <- ggplot(power, aes(x = phi, y = power, color = Model)) +
ylim(0.95, 1) +
geom_line(aes(linetype = Model), size = 1) +
ylab("") +
xlab(expression(phi)) +
ggtitle("Power") +
theme_minimal() +
theme(legend.position = "none") +
scale_linetype_manual(values=c("solid", "longdash", "dotdash"))
# Coverage
cov <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
cov[i, ] <- colMeans(ifelse(mat_lb[[i]] <= beta & mat_ub[[i]] >= beta, 1, 0))
}
cov <- cbind(phi2, cov)
colnames(cov)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
cov <- as.tibble(cov)
cov <- cov %>%
select(-REG) %>%
gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "coverage")
cover <- ggplot(cov, aes(x = phi, y = coverage, color = Model)) +
geom_line(aes(linetype = Model), size = 1) +
ylim(0, 1) +
xlab(expression(phi)) +
ylab("") +
ggtitle("Coverage") +
theme_minimal() +
theme(legend.position = "none") +
scale_linetype_manual(values=c("solid", "longdash", "dotdash"))
# Standard Deviation
sd <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
sd[i, ] <- colSds(mat_coef[[i]])
}
sd <- cbind(phi2, sd)
colnames(sd)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
sd <- as.tibble(sd)
sd_wide <- sd
sd <- sd %>%
select(-REG) %>%
gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "sd")
stddev <- ggplot(sd, aes(x = phi, y = sd, color = Model)) +
geom_line(aes(linetype = Model), size = 1) +
ylim(0, 0.2) +
xlab(expression(phi)) +
ylab("") +
ggtitle("Standard Deviation") +
theme_minimal() +
theme(legend.position = "none") +
scale_linetype_manual(values=c("solid", "longdash", "dotdash"))
# Optimism
mean_se <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
mean_se[i, ] <- colMeans(mat_se[[i]])
}
mean_se <- cbind(phi2, mean_se)
colnames(mean_se)<- c("phi","se_EQ4", "se_LGDV", "se_LGDV2", "se_REG")
mean_se <- as.tibble(mean_se)
opt <- left_join(sd_wide, mean_se, by = "phi")
# opt <- opt %>% mutate(overconf_EQ4 = se_EQ4/EQ4,
#                       overconf_LGDV = se_LGDV/LGDV,
#                       overconf_LGDV2 = se_LGDV2/LGDV2,
#                       overconf_REG = se_REG/REG)
# opt <- opt %>% select(phi, overconf_EQ4, overconf_LGDV, overconf_LGDV2, overconf_REG)
# colnames(opt) <- c("phi","EQ4", "LGDV", "LGDV2", "REG")
# opt <- opt %>%
#   select(-REG) %>%
#   gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "overconf")
opt <- opt %>% mutate(overconf_EQ4 = EQ4/se_EQ4,
overconf_LGDV = LGDV/se_LGDV,
overconf_LGDV2 = LGDV2/se_LGDV2,
overconf_REG = REG/se_REG)
opt <- opt %>% select(phi, overconf_EQ4, overconf_LGDV, overconf_LGDV2, overconf_REG)
colnames(opt) <- c("phi","EQ4", "LGDV", "LGDV2", "REG")
opt <- opt %>%
select(-REG) %>%
gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "overconf")
optimism <- ggplot(opt, aes(x = phi, y = overconf, color = Model)) +
geom_line(aes(linetype = Model), size = 1) +
ylim(0, 2) +
xlab(expression(phi)) +
ylab("") +
ggtitle("Overconfidence") +
theme_minimal() +
theme(legend.position = "none") +
scale_linetype_manual(values=c("solid", "longdash", "dotdash"))
# Combine plots
plot_grid(rootmean, cover, pow, percentbias, stddev,  optimism, nrow = 3)
ggsave("WilkinsRep.pdf", width = 8.5, height = 11, units = "in")
##### ---- Figure 1c and 1d Replication ---- #####
rm(list = ls())
# Load required libraries
library(tidyverse)
library(cowplot)
library(robustbase)
library(matrixStats)
sessionInfo()\
sessionInfo()
